fix: seed FAC supersonic-area-ratio Newton iteration from throat gamma_s - #139

Summary

RocketSolver_solve_supar and RocketSolver_solve_supar_frozen (source/rocket.f90:896,990) hardcode soln%eq_partials(2)%gamma_s as the throat-station seed for the FAC/IAC supersonic-area-ratio Newton iteration. For IAC problems throat_idx=2 so this happened to be correct, but for FAC problems throat_idx=4 (stations: 1=injector, 2=infinity, 3=combustor end, 4=throat), so the infinity station's gamma_s was used instead of the throat's.

Changes

Replace the hardcoded soln%eq_partials(2)%gamma_s with soln%eq_partials(soln%throat_idx)%gamma_s at both sites, matching the pattern already used elsewhere in the file (rocket.f90:661,724) and matching RP-1311 Eq. 6.22, which explicitly states "the value of γₛ is that determined for throat conditions."

Testing

• ctest --test-dir build-dev --output-on-failure — all 15 tests pass, including the FAC-mode rocket_test_suite.test_rocket_fac_only_area_ratios.
• python test/main_interface/test_main.py — 14/14 RP-1311 book-example comparisons pass against published reference values (rtol=1e-2), including example9/example10 (FAC, supar=25,50,75), which exercise exactly this code path.
• Iteration-count check (temporary instrumentation, not included in this diff): Newton-iteration counts for example9/example10 are unchanged before/after the fix (2, 2, 3 iterations) — for these particular cases the infinity- and throat-station gamma_s happen to be close enough that the seed choice doesn't change convergence speed.
• Follow-up sweep, same temporary-instrumentation approach, done in a separate review pass: repeated the iteration-count comparison across 6 additional FAC cases deliberately chosen to maximize the infinity/throat gamma_s gap — H2/O2 at O/F 3–20, RP-1/O2, chamber pressures 10–68 atm, and supar from 2.01 to 150 (largest gap found: ~1.03% relative, at O/F=3, pc=68 atm, supar=100, gamma2=1.23075 vs gammaT=1.24338). Pre-fix and post-fix iteration counts and converged results were identical in every case tried. So while this is a genuine correctness fix — it replaces a physically wrong index with the one RP-1311 Eq. 6.22 specifies — no input found so far (across 8 cases total) demonstrates it changing convergence speed or outcome. The "could hit the 10-iteration cap unconverged" motivation is a theoretical concern based on the code being wrong, not something reproduced numerically; flagging that distinction here rather than overstating the empirical case for this being a robustness fix.

Compatibility / Numerical behavior

• Expected changes (explain and provide validation)

Only the initial-guess seed for a Newton iteration changes; the converged fixed point is unaffected for a converged case (see benchmark results above). Only FAC problems using supar >= 2.0 are affected — IAC problems (throat_idx=2) see no change at all.

Drafted with Claude's assistance

• The station-index mismatch and the fix were confirmed by direct read of rocket.f90, cross-checked against RP-1311 (NASA-CR-19950013764) Eq. 6.21/6.22, which was fetched and read directly (see screenshot in issue history).
• Numerical equivalence was verified by rerunning the full ctest suite and the RP-1311 benchmark comparison script (test/main_interface/test_main.py) against published reference values, not just re-stated from the issue.
• The "no change in iteration count" claim was independently measured by temporarily instrumenting both loop sites, rebuilding, and diffing iteration counts against the pre-fix code — not assumed.
• The follow-up 6-case sweep was added after a request to review whether this PR needed additional test coverage. Rather than write a black-box regression test, the sweep was run first to check whether one was even possible: it confirmed (by actually reverting the fix, rebuilding, and re-running) that no case tried produces an observable difference, so a pass/fail test built on final output couldn't have distinguished the bug from the fix. That finding is reported here instead of a synthetic test.

djkees#10